home *** CD-ROM | disk | FTP | other *** search
/ BCI NET / BCI NET Dec 94.iso / archives / applications / misc / pcal-4.5.lha / pcal / moonphas.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-10-16  |  18.8 KB  |  704 lines
  1. /*
  2.  *  moonphas.c - encapsulates routines used by Pcal for moon phase calculation
  3.  *
  4.  *  The following suite of routines are to calculate the phase of the moon
  5.  *  for a given month, day, year.  They compute the phase of the moon for
  6.  *  noon (UT) on the day requested, the start of the Julian day.
  7.  *
  8.  *  Contents:
  9.  *
  10.  *        calc_phase
  11.  *        find_moonfile
  12.  *        find_phase
  13.  *        read_moonfile
  14.  *
  15.  *  Revision history:
  16.  *
  17.  *    4.5    AWR    11/06/93    accept "opt -[AE]" in the moon file
  18.  *
  19.  *            06/25/93    revise for pre-ANSI math libraries
  20.  *                    without fmod()
  21.  *
  22.  *            04/28/93    restructure function definitions so
  23.  *                    function name appears in first column
  24.  *                    (to facilitate searching for definition
  25.  *                    by name)
  26.  *
  27.  *            04/22/92    eliminated some unused variables and
  28.  *                    calculations
  29.  *
  30.  *    4.4    AWR    01/20/92    use alternate timezone spec (-z)
  31.  *
  32.  *            12/16/91    Revise find_moonfile() for efficiency
  33.  *
  34.  *    4.3    AWR    12/05/91    Search for moon file in Pcal's path
  35.  *                    as well as in calendar file's path
  36.  *
  37.  *            10/25/91    Many changes for support of moon
  38.  *                    phase wildcards and -Z flag
  39.  *
  40.  *    4.11    AWR    08/23/91    Revise is_quarter() to eliminate
  41.  *                    occasional missing or duplicate
  42.  *                    quarter-moons in "-m" mode; add
  43.  *                    gen_phases()
  44.  *
  45.  *    4.1    AWR    08/02/91    Fix bug in julday() where result is
  46.  *                    calculated incorrectly if floor() has
  47.  *                    no prototype
  48.  *
  49.  *      4.01    RLD     03/19/91        Upgraded calc_phase() to accurately
  50.  *                                      calculate the lunar phase for any
  51.  *                                      day without needing to resort to
  52.  *                                      a moon file.
  53.  *
  54.  *    4.0    AWR    03/07/91    Add find_moonfile()
  55.  *
  56.  *            01/15/91    Author: translated PostScript
  57.  *                    routines to C and added moon
  58.  *                    file routines
  59.  */
  60.  
  61. /*
  62.  * Standard headers:
  63.  */
  64.  
  65. #include <stdio.h>
  66. #include <string.h>
  67. #include <ctype.h>
  68. #include <math.h>
  69.  
  70. /*
  71.  * Pcal-specific definitions:
  72.  */
  73.  
  74. #include "pcaldefs.h"
  75. #include "pcalglob.h"
  76. #include "pcallang.h"
  77.  
  78. /*
  79.  * Macros:
  80.  */
  81.  
  82. /*  Astronomical constants. */
  83.  
  84. #define epoch        2444238.5       /* 1980 January 0.0 */
  85.  
  86. /*  Constants defining the Sun's apparent orbit. */
  87.  
  88. #define elonge        278.833540       /* ecliptic longitude of the Sun
  89.                         at epoch 1980.0 */
  90. #define elongp        282.596403       /* ecliptic longitude of the Sun at
  91.                         perigee */
  92. #define eccent      0.016718       /* eccentricity of Earth's orbit */
  93.  
  94. /*  Elements of the Moon's orbit, epoch 1980.0. */
  95.  
  96. #define mmlong      64.975464      /* moon's mean lonigitude at the epoch */
  97. #define mmlongp     349.383063       /* mean longitude of the perigee at the
  98.                                         epoch */
  99. #define mlnode        151.950429       /* mean longitude of the node at the
  100.                         epoch */
  101. #define synmonth    29.53058868    /* synodic month (new Moon to new Moon) */
  102.  
  103. #define PI 3.14159265358979323846  /* assume not near black hole nor in
  104.                         Tennessee ;) */
  105.  
  106.  
  107. /*  Handy mathematical functions. */
  108.  
  109. #define sgn(x) (((x) < 0) ? -1 : ((x) > 0 ? 1 : 0))      /* extract sign */
  110. #ifndef abs
  111. #define abs(x) ((x) < 0 ? (-(x)) : (x))           /* absolute val */
  112. #endif
  113. #define fixangle(a) ((a) - 360.0 * (floor((a) / 360.0)))  /* fix angle      */
  114. #define torad(d) ((d) * (PI / 180.0))              /* deg->rad      */
  115. #define todeg(d) ((d) * (180.0 / PI))              /* rad->deg      */
  116. #define FNITG(x) (sgn (x) * floor (abs (x)))
  117.  
  118. /* fmod() may be absent from pre-ANSI math libraries */
  119.  
  120. #ifdef __STDC__
  121. #define FMOD(x, y)    fmod(x, y)
  122. #else
  123. #define FMOD(x, y)    ((x) - ((long) ((x)/(y)) * (y)))
  124. #endif
  125.  
  126. /* Misc. stuff, to disappear someday along with moon file routines */
  127.  
  128. #define HOUR        12        /* hour of day when phase calculated */
  129. #define PERIOD        synmonth    /* for backward compatibility */
  130.  
  131. /* interpolate phase for day "d" from moon_info array elements "n1" and "n2" */
  132. #define CALC_PHASE(d, n1, n2) \
  133.      moon_info[n1].phase + ((d) - moon_info[n1].doy) * \
  134.     ((moon_info[n2].phase - moon_info[n1].phase) / \
  135.      (moon_info[n2].doy - moon_info[n1].doy))
  136.  
  137. /* generate error message, close file, and quit */
  138. #define ERR_EXIT(msg)    \
  139.     if (1) { ERR(msg); fclose(fp); return FALSE; } else
  140.  
  141. /* day and phase sequence error conditions - cf. read_moonfile() */
  142. #define DAY_TOO_SOON    (nrec > 1 && doy < prevdoy + 6)
  143. #define DAY_TOO_LATE    (doy > prevdoy + 9)
  144. #define WRONG_PHASE    (nrec > 1 && ph != (prevph + 1) % 4)
  145.  
  146.  
  147. /*
  148.  * Globals:
  149.  */
  150.  
  151. typedef struct {
  152.     int    doy;     /* day of year (1..366) */
  153.     int    quarter; /* quarter (MOON_NM, MOON_1Q, etc.) */
  154.     double    phase;     /* moon phase (cycles since new moon prior to 1/1) */
  155. } MOON_INFO;
  156.  
  157. static MOON_INFO moon_info[60];        /* quarter moons for year + dummies */
  158.  
  159.  
  160. /*
  161.  *  Routines to accurately calculate the phase of the moon
  162.  *
  163.  *  Originally adapted from "moontool.c" by John Walker, Release 2.0.
  164.  *
  165.  *      This routine (calc_phase) and its support routines were adapted from
  166.  *      phase.c (v 1.2 88/08/26 22:29:42 jef) in the program "xphoon"        
  167.  *      (v 1.9 88/08/26 22:29:47 jef) by Jef Poskanzer and Craig Leres.
  168.  *      The necessary notice follows...
  169.  *
  170.  *      Copyright (C) 1988 by Jef Poskanzer and Craig Leres.
  171.  *
  172.  *      Permission to use, copy, modify, and distribute this software and its
  173.  *      documentation for any purpose and without fee is hereby granted,
  174.  *      provided that the above copyright notice appear in all copies and that
  175.  *      both that copyright notice and this permission notice appear in
  176.  *      supporting documentation.  This software is provided "as is" without
  177.  *      express or implied warranty.
  178.  *
  179.  *      These were added to "pcal" by RLD on 19-MAR-1991
  180.  */
  181.  
  182.  
  183. /*
  184.  *  julday -- calculate the julian date from input month, day, year
  185.  *      N.B. - The Julian date is computed for noon UT.
  186.  *
  187.  *      Adopted from Peter Duffett-Smith's book `Astronomy With Your
  188.  *      Personal Computer' by Rick Dyson 18-MAR-1991
  189.  */
  190. static double
  191. #ifdef PROTOS
  192. julday(int month,
  193.        int day,
  194.        int year)
  195. #else
  196. julday(month, day, year)
  197.     int month, day, year;
  198. #endif
  199. {
  200.         int mn1, yr1;
  201.         double a, b, c, d, djd;
  202.  
  203.         mn1 = month;
  204.         yr1 = year;
  205.         if ( yr1 < 0 ) yr1 = yr1 + 1;
  206.         if ( month < 3 )
  207.             {
  208.                 mn1 = month + 12;
  209.                 yr1 = yr1 - 1;
  210.             }
  211.         if (( year < 1582 ) ||
  212.             ( year == 1582  && month < 10 ) ||
  213.                 ( year == 1582  && month == 10 && day < 15.0 ))
  214.                     b = 0;
  215.                 else
  216.                 {
  217.                     a = floor (yr1 / 100.0);
  218.                     b = 2 - a + floor (a / 4);
  219.                 }
  220.         if ( yr1 >= 0 )
  221.             c = floor (365.25 * yr1) - 694025;
  222.         else
  223.             c = FNITG ((365.25 * yr1) - 0.75) - 694025;
  224.         d = floor (30.6001 * (mn1 + 1));
  225.         djd = b + c + d + day + 2415020.0;
  226.         return djd;
  227. }
  228.  
  229.  
  230. /*
  231.  *  kepler - solve the equation of Kepler
  232.  */
  233. static double
  234. #ifdef PROTOS
  235. kepler(double m,
  236.        double ecc)
  237. #else
  238. kepler(m, ecc)
  239.     double m, ecc;
  240. #endif
  241. {
  242.     double e, delta;
  243. #define EPSILON 1E-6
  244.  
  245.     e = m = torad(m);
  246.     do {
  247.        delta = e - ecc * sin(e) - m;
  248.        e -= delta / (1 - ecc * cos(e));
  249.     } while (abs(delta) > EPSILON);
  250.     return e;
  251. }
  252.  
  253.  
  254. /*
  255.  *  calc_phase - calculate phase of moon as a fraction:
  256.  *
  257.  *    The argument is the time for which the phase is requested,
  258.  *    expressed as the month, day and year.  It returns the phase of
  259.  *    the moon (0.0 -> 0.99) with the ordering as New Moon, First Quarter,
  260.  *      Full Moon, and Last Quarter.
  261.  *
  262.  *    Converted from the subroutine phase.c used by "xphoon.c" (see
  263.  *      above disclaimer) into calc_phase() for use in "moonphas.c"
  264.  *    by Rick Dyson 18-MAR-1991
  265.  */
  266. double
  267. #ifdef PROTOS
  268. calc_phase(int month,
  269.        int inday,
  270.        int year)
  271. #else
  272. calc_phase(month, inday, year)
  273.     int month, inday, year;
  274. #endif
  275. {
  276.     double Day, N, M, Ec, Lambdasun, ml, MM, Ev, Ae, A3, MmP,
  277.            mEc, A4, lP, V, lPP, MoonAge, pdate;
  278.     static int first = TRUE;
  279.  
  280.     /* get offset from UTC first time through; normalize to +- 1/2 day */
  281.     if (first) {
  282.         utc_offset = FMOD(atof(time_zone), 24.0) / 24.0;
  283.         if (utc_offset < 0.0)
  284.             utc_offset++;
  285.         if (utc_offset > 0.5)
  286.             utc_offset -= 0.5;
  287.         first = FALSE;
  288.         }
  289.  
  290.     /*  need to convert month, day, year into a Julian pdate */
  291.     pdate = julday (month, inday, year) + utc_offset;
  292.  
  293.         /*  Calculation of the Sun's position. */
  294.  
  295.     Day = pdate - epoch;            /* date within epoch */
  296.     N = fixangle((360 / 365.2422) * Day);    /* mean anomaly of the Sun */
  297.     M = fixangle(N + elonge - elongp);      /* convert from perigee
  298.                              co-ordinates to epoch 1980.0 */
  299.     Ec = kepler(M, eccent);            /* solve equation of Kepler */
  300.     Ec = sqrt((1 + eccent) / (1 - eccent)) * tan(Ec / 2);
  301.     Ec = 2 * todeg(atan(Ec));        /* true anomaly */
  302.         Lambdasun = fixangle(Ec + elongp);    /* Sun's geocentric ecliptic
  303.                              longitude */
  304.  
  305.         /*  Calculation of the Moon's position. */
  306.  
  307.         /*  Moon's mean longitude. */
  308.     ml = fixangle(13.1763966 * Day + mmlong);
  309.  
  310.         /*  Moon's mean anomaly. */
  311.     MM = fixangle(ml - 0.1114041 * Day - mmlongp);
  312.  
  313.     /*  Evection. */
  314.     Ev = 1.2739 * sin(torad(2 * (ml - Lambdasun) - MM));
  315.  
  316.     /*  Annual equation. */
  317.     Ae = 0.1858 * sin(torad(M));
  318.  
  319.     /*  Correction term. */
  320.     A3 = 0.37 * sin(torad(M));
  321.  
  322.     /*  Corrected anomaly. */
  323.     MmP = MM + Ev - Ae - A3;
  324.  
  325.     /*  Correction for the equation of the centre. */
  326.     mEc = 6.2886 * sin(torad(MmP));
  327.  
  328.     /*  Another correction term. */
  329.     A4 = 0.214 * sin(torad(2 * MmP));
  330.  
  331.     /*  Corrected longitude. */
  332.     lP = ml + Ev + mEc - Ae + A4;
  333.  
  334.     /*  Variation. */
  335.     V = 0.6583 * sin(torad(2 * (lP - Lambdasun)));
  336.  
  337.     /*  True longitude. */
  338.     lPP = lP + V;
  339.  
  340.     /*  Calculation of the phase of the Moon. */
  341.  
  342.     /*  Age of the Moon in degrees. */
  343.     MoonAge = lPP - Lambdasun;
  344.  
  345.         return (fixangle (MoonAge) / 360.0);
  346. }
  347.  
  348.  
  349. /*
  350.  *  is_quarter - if current phase of moon coincides with quarter moon, return
  351.  *  MOON_NM, MOON_1Q, etc.; otherwise return MOON_OTHER;
  352.  *  
  353.  */
  354. static int
  355. #ifdef PROTOS
  356. is_quarter(double prev,
  357.        double curr,
  358.        double next)
  359. #else
  360. is_quarter(prev, curr, next)
  361.     double prev, curr, next;
  362. #endif
  363. {
  364.     int quarter;
  365.     double phase, diff;
  366.  
  367.     if (curr < prev)    /* adjust phases for 1 -> 0 wraparound */
  368.         curr++;
  369.     if (next < prev)
  370.         next++;
  371.  
  372.     /* if a quarter moon has occurred between "prev" and "next", return
  373.      * TRUE if "curr" is closer to it than "prev" or "next" is    
  374.      */
  375.     for (quarter = 1; quarter <= 4; quarter++) 
  376.         if (prev < (phase = quarter/4.0) && next > phase &&
  377.             (diff = abs(curr - phase)) < phase - prev &&
  378.             diff < next - phase)
  379.             return quarter % 4;        /* MOON_NM == 0 */
  380.  
  381.     return MOON_OTHER;
  382. }
  383.  
  384.  
  385. /*
  386.  * Routines to read moon file and calculate moon phase from data within
  387.  */
  388.  
  389.  
  390. /*
  391.  * make_moonfile - create the base name for the moon file from the supplied
  392.  * template and year; fill it in and return pointer to it
  393.  */
  394. static char *
  395. #ifdef PROTOS
  396. make_moonfile(char *filename,
  397.           char *template,
  398.           int year)
  399. #else
  400. make_moonfile(filename, template, year)
  401.     char *filename;        /* base file name (output) */
  402.     char *template;        /* file name template */
  403.     int year;        /* year to substitute for "%y" in template */
  404. #endif
  405. {
  406.     char *p;
  407.  
  408.     /* copy the template to the file name; replace "%y" (if present)
  409.      * with last two digits of year (as per other year substitutions)
  410.      */
  411.     strcpy(filename, template);
  412.     if ((p = strchr(filename, '%')) && p[1] == 'y') {
  413.         *p++ = '0' + (year / 10) % 10;
  414.         *p   = '0' + year % 10;
  415.     }
  416.  
  417.     return filename;
  418. }
  419.  
  420.  
  421. /*
  422.  * find_moonfile - look for moon file for specified year.  If it exists
  423.  * and is readable, return its full path name; else return NULL.  (There
  424.  * are admittedly ways to do this without attempting to open the file -
  425.  * cf. find_executable() in pcalutil.c - but they may not be portable.)
  426.  */
  427. char *
  428. #ifdef PROTOS
  429. find_moonfile(int year)
  430. #else
  431. find_moonfile(year)
  432.     int year;
  433. #endif
  434. {
  435.     static char filename[STRSIZ];
  436.     char path[STRSIZ], *pathlist[3], moonfile[STRSIZ];
  437.     FILE *fp;
  438.     static int sv_year = 0;
  439.     static char *sv_path = NULL;
  440.  
  441.     /* if -z or -e specified, use algorithm even if moon file present */
  442.     if (tz_flag || ! *datefile)
  443.         return NULL;
  444.  
  445.     /* if year hasn't changed, then return the previous path */
  446.     if (year == sv_year)
  447.         return sv_path;
  448.  
  449.     /* create list of paths for alt_fopen() to search */
  450.     pathlist[0] = mk_path(path, datefile);    /* datefile path */
  451.     pathlist[1] = progpath;            /* program path */
  452.     pathlist[2] = NULL;            /* terminate list */
  453.  
  454.     /* attempt to open the moon file */
  455.     fp = alt_fopen(filename, make_moonfile(moonfile, MOONFILE, year),
  456.                pathlist, "r");
  457.  
  458. #ifdef ALT_MOONFILE
  459.     if (!fp)             /* try again with alternate name */
  460.         fp = alt_fopen(filename, make_moonfile(moonfile, ALT_MOONFILE,
  461.                    year), pathlist, "r");
  462. #endif
  463.  
  464.     /* save year and path for next time around */
  465.     sv_year = year;
  466.     sv_path = fp ? (fclose(fp), filename) : NULL;
  467.  
  468.     return sv_path;
  469. }
  470.  
  471.  
  472. /*
  473.  * read_moonfile - looks for moon data file (in same directory as .calendar
  474.  * or where pcal executable lives - cf. find_moonfile()); if found, reads
  475.  * file, fills in moon_info[] and returns TRUE; if not found (or error
  476.  * encountered), returns FALSE
  477.  */
  478. int
  479. #ifdef PROTOS
  480. read_moonfile(int year)
  481. #else
  482. read_moonfile(year)
  483.     int year;
  484. #endif
  485. {
  486.     char *filename;
  487.     int line, nrec, month, day, hh, mm, dummy, mf_date_style;
  488.     int ph, prevph = MOON_OTHER, doy, prevdoy, n, quarter;
  489.     double phase;
  490.     FILE *fp;
  491.     static char *words[MAXWORD];    /* avoid conflicts with globals */
  492.     static char lbuf[LINSIZ];
  493.  
  494.     /* get name of moon file and attempt to open it */
  495.  
  496.     if ((filename = find_moonfile(year)) == NULL ||
  497.         (fp = fopen(filename, "r")) == NULL) {
  498.         if (DEBUG(DEBUG_MOON) || DEBUG(DEBUG_PATHS))
  499.             FPR(stderr, "No moon file for %d\n", year);
  500.         return FALSE;
  501.     }
  502.  
  503.     if (DEBUG(DEBUG_MOON) || DEBUG(DEBUG_PATHS))
  504.         FPR(stderr, "Using moon file %s\n", filename);
  505.  
  506.     /*
  507.      * Moon file entries are of the form <phase> <date> {<time>}; each
  508.      * is converted below to a moon_info[] record (note that only the
  509.      * initial phase of the moon is directly calculated from <phase>;
  510.      * it is subsequently used only for error checking).  Dummy entries
  511.      * in moon_info[] precede and follow the information read from the
  512.      * moon file; these are used for subsequent interpolation of dates
  513.      * before the first / after the last quarter of the year.
  514.      */
  515.  
  516.     prevdoy = 0;
  517.     line = 0;
  518.     mf_date_style = date_style;    /* unless overridden by "opt -[AE]" */
  519.  
  520.     for (nrec = 1; getline(fp, lbuf, &line); nrec++) {
  521.  
  522.         /* special check for "opt -[AE]" line */
  523.         n = loadwords(words, lbuf);
  524.  
  525.         if (n == 2 && date_type(words[0], &dummy, &dummy) == DT_OPT) {
  526.             if (words[1][0] == '-') {
  527.                 char c = words[1][1];
  528.                     if (c == F_USA_DATES || c == F_EUR_DATES) {
  529.                     mf_date_style = c == F_USA_DATES ?
  530.                             USA_DATES : EUR_DATES;
  531.                     nrec--;        /* skip this line */
  532.                     continue;
  533.                 }
  534.             }
  535.         }
  536.  
  537.         /* ensure line is recognizable */
  538.         if (n < 2 || (ph = get_phase(words[0])) == MOON_OTHER)
  539.             ERR_EXIT(E_INV_LINE);
  540.     
  541.         if (nrec == 1)            /* phase at initial quarter */
  542.             quarter = ph == MOON_NM ? 4 : ph;
  543.  
  544.         /* extract the month and day fields (in appropriate order) */
  545.  
  546.         (void) split_date(words[1],
  547.                   mf_date_style == USA_DATES ? &month : &day,
  548.                   mf_date_style == USA_DATES ? &day : &month,
  549.                   NULL);
  550.  
  551.         /* validate the date and phase */
  552.  
  553.         if (!is_valid(month, day, year))    /* date OK? */
  554.             ERR_EXIT(E_INV_DATE);
  555.  
  556.         doy = DAY_OF_YEAR(month, day, year);    /* in sequence? */
  557.         if (DAY_TOO_SOON || DAY_TOO_LATE || WRONG_PHASE)
  558.             ERR_EXIT(E_DATE_SEQ);
  559.  
  560.         prevdoy = doy;            /* save for sequence check */
  561.         prevph = ph;
  562.  
  563.         /* calculate moon phase, factoring in time (if present) */
  564.  
  565.         phase = 0.25 * quarter++;
  566.         if (n > 2) {            /* extract hour and minute */
  567.             (void) split_date(words[2], &hh, &mm, NULL);
  568.             phase += (HOUR - (hh + (mm / 60.0))) / (24 * PERIOD); 
  569.         }
  570.         moon_info[nrec].doy = doy;    /* enter day and phase */
  571.         moon_info[nrec].quarter = ph % 4;
  572.         moon_info[nrec].phase = phase;
  573.  
  574.         if (DEBUG(DEBUG_MOON))
  575.             FPR(stderr, "moon_info[%2d]: %3d %d %6.3f\n", nrec,
  576.                 moon_info[nrec].doy, moon_info[nrec].quarter,
  577.                 moon_info[nrec].phase);
  578.     }
  579.  
  580.     /* check to see that file is all there */
  581.  
  582.     doy = YEAR_LEN(year) + 1;    /* day after end of year */
  583.     if (DAY_TOO_LATE)
  584.         ERR_EXIT(E_PREM_EOF);
  585.  
  586.     /* extrapolate dummy entries from nearest lunar month */
  587.  
  588.     moon_info[nrec].doy = doy;    /* day after end of year */
  589.     moon_info[nrec].quarter = MOON_OTHER;
  590.     moon_info[nrec].phase = CALC_PHASE(doy, nrec-5, nrec-1);
  591.     if (DEBUG(DEBUG_MOON))
  592.         FPR(stderr, "moon_info[%2d]: %3d %d %6.3f\n", nrec,
  593.             moon_info[nrec].doy, moon_info[nrec].quarter,
  594.             moon_info[nrec].phase);
  595.  
  596.     moon_info[0].doy = 0;        /* day before start of year */
  597.     moon_info[0].quarter = MOON_OTHER;
  598.     moon_info[0].phase = CALC_PHASE(0, 1, 5);
  599.     if (DEBUG(DEBUG_MOON))
  600.         FPR(stderr, "moon_info[%2d]: %3d %d %6.3f\n", 0,
  601.             moon_info[0].doy, moon_info[0].quarter,
  602.             moon_info[0].phase);
  603.     
  604.     fclose(fp);
  605.     return TRUE;
  606. }
  607.  
  608.  
  609. /*
  610.  * gen_phases - fill array with moon phases for all days in month/year (plus
  611.  * extra entries at beginning and end for last day of previous month and
  612.  * first day of next month, respectively)
  613.  */
  614. static void
  615. #ifdef PROTOS
  616. gen_phases(double phase[],
  617.        int month,
  618.        int year)
  619. #else
  620. gen_phases(phase, month, year)
  621.     double phase[];
  622.     int month, year;
  623. #endif
  624. {
  625.     int day, len;
  626.     DATE date;
  627.  
  628.     /* start with moon phase for last day of previous month */
  629.     MAKE_DATE(date, month, 0, year);
  630.     normalize(&date);
  631.     phase[0] = calc_phase(date.mm, date.dd, date.yy);
  632.  
  633.     /* add the moon phases for all days in the current month */
  634.     for (day = 1, len = LENGTH_OF(month, year); day <= len; day++)
  635.         phase[day] = calc_phase(month, day, year);
  636.  
  637.     /* append the moon phase for the first day of next month */
  638.     MAKE_DATE(date, month, len + 1, year);
  639.     normalize(&date);
  640.     phase[len + 1] = calc_phase(date.mm, date.dd, date.yy);
  641. }
  642.  
  643.  
  644. /*
  645.  * find_phase - look up phase of moon in moon phase file (if possible);
  646.  * otherwise calculate it using calc_phase() above.  Sets *pquarter to
  647.  * MOON_NM, MOON_1Q, etc. if quarter moon, MOON_OTHER if not
  648.  */
  649. double
  650. #ifdef PROTOS
  651. find_phase(int month,
  652.        int day,
  653.        int year,
  654.        int *pquarter)
  655. #else
  656. find_phase(month, day, year, pquarter)
  657.     int month, day, year;
  658.     int *pquarter;
  659. #endif
  660. {
  661.     static int sv_year = 0, sv_month = 0;
  662.     static int use_file;
  663.     static double mphase[33];    /* 31 days + 2 dummies */
  664.     int i, doy;
  665.     double phase;
  666.  
  667.     if (year != sv_year) {        /* look for file for new year */
  668.         use_file = read_moonfile(year);
  669.     }
  670.  
  671.     if (! use_file) {        /* no file - calculate phase */
  672.         /* new month?  fill mphase[] with moon phases */
  673.         if (month != sv_month || year != sv_year) {
  674.             gen_phases(mphase, month, year);
  675.             sv_month = month;
  676.             sv_year = year;
  677.         }
  678.  
  679.         phase = mphase[day];
  680.         *pquarter = is_quarter(mphase[day-1], phase, mphase[day+1]);
  681.         return phase;
  682.     }
  683.  
  684.     /* moon file found - use the data read by read_moonfile() */
  685.  
  686.     sv_year = year;                /* avoid re-reading same file */
  687.     doy = DAY_OF_YEAR(month, day, year);
  688.     for (i = 1; doy > moon_info[i].doy; i++)    /* find interval */
  689.         ;
  690.  
  691.     /* if day appears in table, return exact value; else interpolate */
  692.  
  693.     if (doy == moon_info[i].doy) {
  694.         *pquarter = moon_info[i].quarter;
  695.         phase = moon_info[i].phase;
  696.     } else {
  697.         *pquarter = MOON_OTHER;
  698.         phase = CALC_PHASE(doy, i-1, i);
  699.     }
  700.  
  701.     return phase - (int)phase;            /* 0.0 <= phase < 1.0 */
  702. }
  703.  
  704.